{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Multivariable pairing (RGA)\n", "\n", "For a 2$\\times$2 system, we have 2 choices of pairing variables for distributed control:\n", "\n", "\n", "\n", "\n", "\n", "
DiagonalOff-diagonal
$$G_{cd} = \\left[\\begin{array}{cc} G_{c1} & 0 \\\\ 0 & G_{c2} \\end{array}\\right]$$$$G_{co} = \\left[\\begin{array}{cc} 0 & G_{c2} \\\\ G_{c1} & 0 \\end{array}\\right]$$
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bristol developed the Relative Gain Array to determine good pairings based on only the plant transfer function matrix $G_p$. The elements of the RGA are defined as\n", "\n", "$$\\lambda_{ij} \\triangleq \\frac{(\\partial y_i/\\partial u_j)_u}{(\\partial y_i/\\partial u_j)_y}= \\frac{\\text{open loop gain}}{\\text{closed loop gain}}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could build $\\Lambda$ by direct evaluation of the above derivatives near some point given a time-domain model, but if we already have a transfer function model, we can evaluate the steady-state gain matrix $K$ by using the final value theorem." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import sympy\n", "sympy.init_printing()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "s = sympy.Symbol('s')" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def fopdt(k, theta, tau):\n", " return k*sympy.exp(-theta*s)/(tau*s + 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the system from example 16.5" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\\left[\\begin{matrix}- \\frac{2 e^{- s}}{10 s + 1} & \\frac{1.5 e^{- s}}{s + 1}\\\\\\frac{1.5 e^{- s}}{s + 1} & \\frac{2 e^{- s}}{10 s + 1}\\end{matrix}\\right]$$" ], "text/plain": [ "⎡ -s -s ⎤\n", "⎢-2⋅ℯ 1.5⋅ℯ ⎥\n", "⎢──────── ─────── ⎥\n", "⎢10⋅s + 1 s + 1 ⎥\n", "⎢ ⎥\n", "⎢ -s -s ⎥\n", "⎢1.5⋅ℯ 2⋅ℯ ⎥\n", "⎢─────── ────────⎥\n", "⎣ s + 1 10⋅s + 1⎦" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "G_p = sympy.Matrix([[fopdt(-2, 1, 10), fopdt(1.5, 1, 1)],\n", " [fopdt(1.5, 1, 1), fopdt(2, 1, 10)]])\n", "G_p" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Unfortunately sympy cannot calculate limits on matrix expressions" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "#K = sympy.limit(G_p, s, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But we can apply a function to the elements:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def gain(G):\n", " return sympy.limit(G, s, 0)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "K = G_p.applyfunc(gain)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\\left[\\begin{matrix}-2 & 1.5\\\\1.5 & 2\\end{matrix}\\right]$$" ], "text/plain": [ "⎡-2 1.5⎤\n", "⎢ ⎥\n", "⎣1.5 2 ⎦" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "K" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then calculate $\\Lambda = K \\otimes H$ where $H = (K^{-1})^{T}$:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$\\left[\\begin{matrix}0.64 & 0.36\\\\0.36 & 0.64\\end{matrix}\\right]$$" ], "text/plain": [ "⎡0.64 0.36⎤\n", "⎢ ⎥\n", "⎣0.36 0.64⎦" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Lambda = K.multiply_elementwise(K.inv().transpose())\n", "Lambda" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can do the same calculation (faster) using numpy:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "import numpy" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def fopdt(k, theta, tau):\n", " return k*numpy.exp(-theta*s)/(tau*s + 1)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "s = 0" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "K = numpy.matrix([[fopdt(-2, 1, 10), fopdt(1.5, 1, 1)],\n", " [fopdt(1.5, 1, 1), fopdt(2, 1, 10)]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `.A` attribute in matrices is the matrix as a `numpy.array`, which multiplies elementwise by default." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.64, 0.36],\n", " [0.36, 0.64]])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "K.A*K.I.T.A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The numpy developers recommended that you should use `numpy.array` instead of `numpy.matrix` as much as possible. I find this makes the notation harder to read:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "K = numpy.array([[fopdt(-2, 1, 10), fopdt(1.5, 1, 1)],\n", " [fopdt(1.5, 1, 1), fopdt(2, 1, 10)]])" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0.64, 0.36],\n", " [0.36, 0.64]])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "K*numpy.linalg.inv(K).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulation results\n", "\n", "Let's simulate this system to get an idea of how the control works out" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "import tbcontrol\n", "tbcontrol.expectversion('0.1.4')\n", "from tbcontrol import blocksim" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "import numpy" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "N = 2" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "G = {}" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "gains = [[-2, 1.5], \n", " [1.5, 2]]\n", "taus = [[10, 1], \n", " [1, 10]]\n", "delays = [[1, 1], \n", " [1, 1]]" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "for inp in range(N):\n", " for outp in range(N):\n", " G[(outp, inp)] = blocksim.LTI(f\"G_{outp}_{inp}\", f\"u_{inp}\", f\"yp_{inp}_{outp}\", \n", " gains[outp][inp], [1, taus[outp][inp]], delays[outp][inp])" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "inputs = {'ysp_0': blocksim.step(),\n", " 'ysp_1': blocksim.step(starttime=50)}" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "sums = {f'y_{outp}': [f\"+yp_{inp}_{outp}\" for inp in range(N)] for outp in range(N)}\n", "for i in range(N):\n", " sums[f'e_{i}'] = [f'+ysp_{i}', f'-y_{i}']" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'y_0': ['+yp_0_0', '+yp_1_0'],\n", " 'y_1': ['+yp_0_1', '+yp_1_1'],\n", " 'e_0': ['+ysp_0', '-y_0'],\n", " 'e_1': ['+ysp_1', '-y_1']}" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sums" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "def simulate(auto1=True, K1=-1, tauI1=10, auto2=True, K2=0.5, tauI2=10):\n", " controllers = {'Gc_0': blocksim.PI('Gc_0', 'e_0', 'u_0', K1, tauI1),\n", " 'Gc_1': blocksim.PI('Gc_1', 'e_1', 'u_1', K2, tauI2)}\n", "\n", " controllers['Gc_0'].automatic = auto1\n", " controllers['Gc_1'].automatic = auto2\n", "\n", " ts = numpy.arange(0, 100, 0.125)\n", "\n", " diagram = blocksim.Diagram(list(G.values()) + list(controllers.values()), sums, inputs)\n", "\n", " result = diagram.simulate(ts)\n", "\n", "# plt.figure()\n", "# plt.plot(ts, result['u_0'])\n", "# plt.plot(ts, result['u_1'])\n", "\n", " plt.figure()\n", " plt.plot(ts, result['y_0'])\n", " plt.plot(ts, result['ysp_0'])\n", " plt.plot(ts, result['y_1'])\n", " plt.plot(ts, result['ysp_1'])" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "from ipywidgets import interact" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "f544acf5cc3f463091c80c5472d35b15", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(Dropdown(description='auto1', options=(True, False), value=True), FloatSlider(value=-1.0…" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "interact(simulate, \n", " auto1=[True, False], K1=(-2., 0), tauI1=(1., 50), \n", " auto2=[True, False], K2=(0., 2), tauI2=(1., 50))" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }